First principles simulations of direct coexistence of solid and liquid aluminium 
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First principles calculations based on density functional theory, with generalised gradient corrections 
and ultrasoft pseudopotentials, have been used to simulate solid and liquid aluminium in direct 
coexistence at zero pressure. Simulations have been carried out on systems containing up to 1000 
atoms for 15 ps. The points on the melting curve extracted from these simulations are in very good 
agreement with previous calculations, which employed the same electronic structure method but 
used an approach based on the explicit calculation of free energies [1]. 
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The calculation of melting properties of materials us- 
ing computer simulations has a long history. Two main 
approaches have traditionally been used. The first is 
based on the direct simulation of solid and liquid in co- 
existence [2-5,17]. The second on the calculation of the 
free energy of sohd and hquid [7-10,18], with the melting 
point (p, T) determined by the condition of equality of 
the Gibbs free energies of hquid and solid, Giiq(p, T) = 
Gso\ (p, T) . The two approaches must clearly give the 
same answer once all the sources of errors have been 
brought under control. These works were originally car- 
ried out using classical model potentials, with parameters 
usually adjusted to reproduce some known experimental 
properties of the material, or fitted to ab-initio calcula- 
tions. The advantage of using these classical potentials is 
that they are relatively simple, so that computer simula- 
tions could easily be carried out on large systems and for 
long time. The main disadvantage, however, is that the 
transferability of these potentials is not always guaran- 
teed, so that the accuracy of the predictions is sometimes 
questionable. In 1995 Sugino and Car [11] showed that 
it was actually possible to use density functional theory 
(DFT) [12] techniques to calculate the melting temper- 
ature of materials truly from first principles, i.e without 
relying on any adjustable parameters. They chose the 
melting of silicon as a test case, and using the local den- 
sity approximation (LDA) they calculated the Gibbs free 
energy of solid and liquid using thermodynamic integra- 
tion. This is a well known statistical mechanics tech- 
nique to compute free energy differences between two sys- 
tems [13]. The idea is to use a simple system for which 
the free energy is known as a reference system, and then 
compute the free energy difference between the ab-initio 
and the reference system, which is equal to the reversible 
work done in adiabatically switching the potential energy 
from one system to the other. An attractive feature of 
this method is that the final result is totally independent 
on the choice of the reference system. In practice, how- 
ever, this method can only work if one is able to find a 
reference system which is as close as possible to the ah 
initio system, so that the computational effort needed 
to calculate the free energy difference between the two 



systems is reduced to the minimum. Using this method, 
Sugino and Car found a zero pressure melting tempera- 
ture about 20% lower than the experimental datum. This 
may seem not particularly good, but it has to be remem- 
bered that no adjustable parameters were inputed in the 
calculations, nor any experimental data other than the 
Plank's constant and the mass and charge of the elec- 
tron. Recently, we have argued [14] that this non-perfect 
agreement between LDA results and experiments is prob- 
ably due to non-cancelling LDA errors between the solid 
and liquid. The reason is that the two phases have very 
different properties: the liquid is six-fold coordinated and 
the solid four-fold coordinated, the liquid is a metal and 
the solid an insulator, so it is not surprising that the 
LDA errors in the two phases may be different. We have 
repeated the LDA calculations of Sugino and Car using 
similar techniques and found the same melting tempera- 
ture. However, we have also shown that if the generalised 
gradient corrections (GGA) are used, then the melting 
temperature of silicon comes out in much closer agree- 
ment with the experiments [14]. 

Since the work of Sugino and Car [11] a number of 
first principles based calculations of melting points and 
melting curves have followed. In some cases the free en- 
ergy approach has been used [15,16,19,20,1], while other 
methods relied on fitting a model potential to first princi- 
ples simulations and calculating melting properties with 
the model potential [3-5] . The advantage of the free en- 
ergy approach is that it is unbiased, provided all sources 
of technical errors are brought under control. A disad- 
vantage of the method is that it is intrinsically complex. 
On the other hand, the coexistence approach is relatively 
simple to apply, but has the main disadvantage of rely- 
ing on good quality fitting and, more importantly, trans- 
ferability of the model for at least a simultaneous good 
description of solid and liquid. Recently, we have shown 
that provided that the model potential is reasonably close 
to the ab-initio system [21], it is possible to correct for 
the (small) differences between the model and the full 
ab-initio system using a perturbational approach to ther- 
modynamic integration, and that once the corrections are 
applied the results coincide with those obtained using the 
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free energy approach [22] . 

The coexistence approach has been used extensively in 
the past to calculate the melting properties of various ma- 
terials. In the constant volume constant internal energy 
[NVE ensemble) approach to the method it has been 
shown that liquid and solid can coexist for long time, 
provided V and E arc appropriately chosen [2,22]. The 
average value of the pressure p and temperature T over 
the coexisting period then give a point on the melting 
curve. Size effects have also been studied quite exten- 
sively, and it was shown that correct results can be ob- 
tained in systems containing more than 500 atoms [5,6]. 

Here I have exploited recent advances in computer 
power and algorithms developments to use direct DFT 
calculations for simulating solid and liquid aluminium in 
coexistence near zero pressure. This work combines the 
simplicity of the coexistence approach with the accuracy 
provided by DFT, and in some respect represents a shift 
of paradigm, whereby the main effort is transferred from 
the human to the computer. The temperature at which 
solid and liquid coexist are found to be in good agreement 
with our previous results obtained with the free energy 
approach [1], so that these results also provide additional 
evidence that our techniques to calculate free energies are 
sound. 

The calculations have been performed with the VASP 
code [23] , with the implementation of an efhcient extrap- 
olation of the charge density [24], ultrasoft pseudopo- 
tentials [25] with a plane wave cutoff of 130 eV, and 
generalised gradient corrections. The simulations have 
been performed in the NVE ensemble on systems con- 
taining 1000 atoms (5 x 5 x 10 cubic supercell), with 
a time step of 3 fs and a convergency threshold on the 
total energy of 2 x 10~^ eV/atom. With these prescrip- 
tions the drift in the microcanonical total energy was 
less than 0.3 K/ps. The total length of the simulations 
were typically 15 ps. The volume per atom was cho- 
sen to be V = 18.5 A'\ which is close to the average 
of the volumes of solid and liquid at the zero pressure 
melting point calculated in our previous work [1]. Elec- 
tronic excitations have been included within the frame- 
work of finite temperature DFT. The simulations have 
been performed using the F point only, and spot checked 
with Monkhorst-Pack [26] (2 x 2 x 1) and (4 x 4 x 2) 
k-point grids. Calculations with the F point predict a 
small non-hydrostatic stress tensor with a difference of 
about Pz — {px '^Py)l'2- = —2 kB between the compo- 
nents of the stress parallel to the solid-liquid interface, 
Px and Py, and that perpendicular to the interface, Pz- 
The pressure p= {Px+Py+Pz)/'^ ''^^ essentially exact and 
the three off-diagonal components of the stress tensor 
fluctuate around zero average, so that there is no shear 
stress on the cell. With F the total energy is wrong by 
?» 5 meV/atom, and it is likely that the error is almost 
equally shared by the liquid and the solid parts, so that 
the resulting error on the melting temperature is most 



probably negligible [27]. A correction term of 2.7 kB due 
to the lack of convergency with respect to the plane wave 
cutoff has been added to the calculated pressures. The 
systems have been monitored by inspecting the average 
number density in slices of the cell taken parallel to the 
boundary between solid and liquid. In the solid region 
this number is a periodic function of the slice number, 
and in the liquid part it fluctuates randomly around its 
average value. 

The zero pressure crystal structure of aluminium is 
face-centred-cube (fee), so I have assumed that melting 
occurs from this structure. To prepare the system I have 
used the inverse power classical potential employed in 
Ref. [1]. This model has been tuned to the same ab-initio 
system used here, so it is a good starting point for the 
present calculations. The preparation procedure follows 
Ref. [22] . A perfect crystal is initially thermalized at 800 
K, then the simulation is stopped, half of the atoms are 
clamped and the other half are freely evolved at very 
high temperature until melting occurs, then the liquid is 
thermalized back at 800 K [28] . At this point the system 
is being freely evolved in the NVE ensemble using DFT. 

As explained in Ref. [22] , for each chosen volume there 
is a whole range of internal energies for which different 
amount of solid and liquid are in coexistence. Provided 
that the energy is not too low (high) so that that the 
whole system freezes (melts), a whole piece of melting 
curve can be obtained for any fixed volume. In this work I 
have performed 3 simulations with 1000 atoms, all at the 
same volume but with different amounts of total energy, 
which therefore provide three distinct points on the melt- 
ing curve. In all simulations coexistence was obtained for 
the whole length of the runs. In Fig. 1 I display the den- 
sity profile, calculated by dividing the simulations cell 
into 100 slices parallel to the solid-liquid interface, corre- 
sponding to the last configuration of one of these simula- 
tions. Fig. 2 contains the calculated temperature (upper 
panel) and pressure (lower panel) for this particular sim- 
ulation, and shows that they oscillate stably around their 
average values T « 820 K and p w 5.5 kB [29]. In order 
to test the reliability of the lengths of these simulations I 
have performed additional runs on system containing 512 
atoms (4 X 4 X 10 cubic supercell), which could be simu- 
lated for up to 40 ps. Temperature and pressure from one 
of these runs are displayed in Fig. 3. It is clear that all 
the information needed to extract useful values for p and 
T is contained in any time window of rs 5 — 10 ps, which 
provides confidence that the simulations for the large sys- 
tems are long enough. A more detailed inspection of the 
figures reveals the presence of anti-correlated oscillations 
in pressure and temperature, which correspond to fluc- 
tuations in the total amount of solid and liquid in the 
system. These fluctuations seem absent in the simula- 
tions with 1000 atoms, but they would probably develop 
if the runs could be extended. In any case, they do not 
affect the average value of p and T. 
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It has been pointed out that non-hydrostaic conditions 
artificially raise the free energy of the solid, and there- 
fore lower the melting temperature [6]. To investigate 
the effect on the melting temperature due to the present 
non-hydrostatic conditions I have used the model po- 
tential used in Ref. [1]. The model showed a difference 
Pz — {Px ~^Py)/'^ = 0.5 kB when simulated with the same 
cell as the ab-initio system with 1000 atoms. To estimate 
the systematic error on the melting temperature I have 
performed an additional simulation using a slightly elon- 
gated cell at the same volume, so that non-hydrostaticity 
was reduced to less than 0.1 kB. The decrease of melting 
temperature in this second simulation was less than 10 
K. A similar effect is expected for the ab-initio system, so 
that I estimate a systematic error due to non-hydrostatic 
conditions to be at worse of the order of 20 K. 

In Fig. 4 I report the values of p and T calculated here 
in comparison with the low end of the melting curve cal- 
culated in Ref. [1] The agreement between the present 
finding with 1000 atoms and the free energy approach 
results is extremely good. The results obtained from the 
simulations with 512 atoms display melting temperatures 
higher by about 50 K, although they are still compatible 
with the results obtained with the free energy approach 
within the combined error bars. In order to test if the 
cause of this size effect was due to inadequate k-point 
sampling, I have performed a simulation using a (2 x 2 x 1) 
k-point grid on a system with 512 atoms for « 6 ps. 
The length of this simulation is too short to draw defi- 
nite conclusions, but it indicates that the pressure may 
be underestimated by « 1 — 2 kB and the temperature 
over-predicted by « 20 K. This would bring the results 
in somewhat better agreement with those obtained with 
the larger system, and also with the calculations based 
on the free energy approach (which were fully converged 
with respect to size and k-point sampling) [1]. Notice 
that these results do not agree perfectly with the experi- 
mental zero pressure melting temperature of aluminium 
(933 K) [30]. We have argued in our previous work [1] 
that this disagreement is due to inadequacy of the GGA 
to predict the vibrational properties of the solid. We have 
also suggested that this inaccuracy can be rationalised 
in terms of the errors in the equation of state, and de- 
vised a simple way to correct for it. In Table 1 I report 
the calculated GGA lattice constant, the bulk modulus 
and the cohesive energy compared with the experimental 
data. Following Gaudoin et al. [32], I also report the ex- 
perimental values adjusted for the absence of zero point 
motion effects, which is missing in the present calculated 
values. It is evident that the GGA predicts a zero pres- 
sure lattice constant for Al which is too large, or, which 
means that at the correct lattice constant GGA predicts 
a positive pressure of about 16 kB. It follows that the 
GGA zero pressure melting point is actually the melting 
point at a negative pressure of about -16 kB, which is 
about 120-130 K lower. 



In summary, I have shown here that it has become pos- 
sible to directly simulate solid and liquid in coexistence 
using density functional theory techniques. This method 
combines the advantages of being relative easy to use and 
provide DFT accuracy. It is inevitably computationally 
very intensive, and these calculations have only been pos- 
sible thanks to an access onto the 3.5 Teraflop/s machine 
(IBM p690 Regatta with 1280 processors). The cost of 
each molecular dynamics step for the 1000 atoms system 
was between 3 and 4 minutes on 128 processors of the 
machine. However, as large computer resources become 
routinely available this method should find widespread 
applicability in the future. The present results are in 
very good agreement with our previous findings based 
on the direct calculations of free energies, and therefore 
also support the reliability of those techniques. 
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Experiment 


GGA 


E 


3.39 (3.43) 


3.43 


B 


76 (81) 


73 


ao 


4.049 (4.022) 


4.05 



TABLE 1. Experimental values for the cohesive energy 
E" (eV/atom) [31], the bulk modulus B (GPa) [30] and 
the lattice constant ao (A) [30]. In parenthesis axe the ad- 
justed values for the effect of zero point motion [32]. The 
calculated GGA values have been obtained from a fit to a 
Birch- Murnaghan equation of state [33], and they do not in- 
clude zero point motion effects. 
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FIG. 1. Density profile in a simulation of solid and liquid Al 
coexisting at zero pressure. The system is divided in slices of 
equal thickness (0.42 A) parallel to the solid-liquid interface, 

and the graph shows number of atoms in each slice. Simula- 
tions were performed on a system of 1000 atoms using density 
functional theory. 
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FIG. 2. Time variation of temperature (upper panel) and 
pressure (lower panel) during a simulation of solid and liquid 
Al in coexistence. Simulations were performed on a system of 

1000 atoms with density functional theory. Gray lines: actual 
data; black lines: running averages over a 0.75 ps period. 
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FIG. 3. Time variation of temperature (upper panel) and 
pressure (lower panel) during a simulation of solid and liquid 
Al in coexistence for a system containing 512 atoms. Gray 
lines: actual data; black lines: running averages over a 0.75 
ps period. 
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FIG. 4. Temperatures and pressures at which liquid and 
solid coexist in simulations containing 1000 atoms (circles) 
and 512 atoms (triangles). The solid line is the lower end of 
the melting curve calculated using the free energy approach 
in Ref. [1] (see text), light dashed lines represent error bars. 



5 



